# Importing Libraries

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm, qmc
from tqdm import tqdm
from numba import njit
from joblib import Parallel, delayed
%matplotlib inline

@dataclass
class MonteCarloOptionPricing:
    """Monte Carlo Option Pricing Engine"""
    
    S0: float = field(default=100.0, metadata={"description": "Initial stock price", "gt": 0})
    K: float = field(default=100.0, metadata={"description": "Strike price", "gt": 0})
    r: float = field(default=0.05, metadata={"description": "Risk-free rate", "ge": 0, "le": 1})
    sigma: float = field(default=0.2, metadata={"description": "Volatility", "gt": 0, "le": 1})
    T: float = field(default=1.0, metadata={"description": "Time to expiration (years)", "gt": 0, "le": 10})
    N: int = field(default=10000, metadata={"description": "Number of simulations", "gt": 0, "le": 1_000_000})
    ts: int = field(default=252, metadata={"description": "Number of timesteps", "gt": 0})
    seed: Optional[int] = field(default=None, metadata={"description": "Random seed for reproducibility"})
    
    # Internal attributes
    S: Optional[np.ndarray] = field(init=False, default=None, repr=False)
    discount_factor: float = field(init=False)
    
    def __post_init__(self):
        # Set the discount factor
        self.discount_factor = np.exp(-self.r * self.T)
        # Set the random seed if provided
        if self.seed is not None:
            np.random.seed(self.seed)
    
    def simulate_euler_maruyama(self) -> np.ndarray:
        """Simulate stock price paths using the Euler-Maruyama scheme."""
        dt = self.T / self.ts
        drift = (self.r - 0.5 * self.sigma**2) * dt
        diffusion = self.sigma * np.sqrt(dt)
        Z = np.random.standard_normal((self.N, self.ts))
        log_S = np.log(self.S0) + np.cumsum(drift + diffusion * Z, axis=1)
        S = np.exp(log_S)
        S = np.hstack((np.full((self.N, 1), self.S0), S))
        self.S = S
        return S
    
    def simulate_antithetic_variates(self) -> np.ndarray:
        """Simulate stock price paths using Antithetic Variates for variance reduction."""
        dt = self.T / self.ts
        drift = (self.r - 0.5 * self.sigma**2) * dt
        diffusion = self.sigma * np.sqrt(dt)
        
        M_half = self.N // 2
        Z = np.random.standard_normal((M_half, self.ts))
        Z_antithetic = -Z
        Z_full = np.vstack((Z, Z_antithetic))
        
        log_S = np.log(self.S0) + np.cumsum(drift + diffusion * Z_full, axis=1)
        S = np.exp(log_S)
        S = np.hstack((np.full((self.N, 1), self.S0), S))
        self.S = S
        return S
    
    def simulate_quasi_monte_carlo(self) -> np.ndarray:
        """Simulate stock price paths using Quasi-Monte Carlo with Sobol sequences."""
        sampler = qmc.Sobol(d=self.ts, scramble=True)
        # Sobol sequences require M to be a power of 2; generate more and truncate
        M_power = 2**int(np.ceil(np.log2(self.N)))
        sample = sampler.random(M_power)
        # Inverse transform to standard normal
        Z = norm.ppf(sample)
        Z = Z[:self.N]  # Truncate to N samples
        
        dt = self.T / self.ts
        drift = (self.r - 0.5 * self.sigma**2) * dt
        diffusion = self.sigma * np.sqrt(dt)
        
        log_S = np.log(self.S0) + np.cumsum(drift + diffusion * Z, axis=1)
        S = np.exp(log_S)
        S = np.hstack((np.full((self.N, 1), self.S0), S))
        self.S = S
        return S
    
    def simulate_milstein(self) -> np.ndarray:
        """Simulate stock price paths using Milstein’s method for higher accuracy."""
        dt = self.T / self.ts
        drift = (self.r - 0.5 * self.sigma**2) * dt
        diffusion = self.sigma * np.sqrt(dt)
        
        Z = np.random.standard_normal((self.N, self.ts))
        
        # Euler-Maruyama step
        log_S = np.log(self.S0) + np.cumsum(drift + diffusion * Z, axis=1)
        
        # Milstein correction
        milstein_correction = 0.5 * self.sigma**2 * (Z**2 - 1) * dt
        log_S += np.cumsum(milstein_correction, axis=1)
        
        S = np.exp(log_S)
        S = np.hstack((np.full((self.N, 1), self.S0), S))
        self.S = S
        return S
    
    def simulate_paths(self, method: str = 'Euler-Maruyama') -> np.ndarray:
        """Simulate stock price paths using the specified method."""
        if method == 'Euler-Maruyama':
            return self.simulate_euler_maruyama()
        elif method == 'Antithetic Variates':
            return self.simulate_antithetic_variates()
        elif method == 'Quasi-Monte Carlo':
            return self.simulate_quasi_monte_carlo()
        elif method == 'Milstein':
            return self.simulate_milstein()
        else:
            raise ValueError(f"Unsupported simulation method: {method}")
    
    def european_call_price_bs(self) -> float:
        """Calculate the analytical Black-Scholes price for a European Call option."""
        d1 = (np.log(self.S0 / self.K) + (self.r + 0.5 * self.sigma**2) * self.T) / (self.sigma * np.sqrt(self.T))
        d2 = d1 - self.sigma * np.sqrt(self.T)
        price = self.S0 * norm.cdf(d1) - self.K * np.exp(-self.r*self.T) * norm.cdf(d2)
        return price
    
    def european_call_payoff(self) -> np.ndarray:
        """Calculate payoff for European Call option."""
        return np.maximum(self.S[:, -1] - self.K, 0)
    
    def asian_call_payoff(self) -> np.ndarray:
        """Calculate payoff for Asian Call option."""
        S_avg = self.S[:, 1:].mean(axis=1)
        return np.maximum(S_avg - self.K, 0)
    
    def asian_put_payoff(self) -> np.ndarray:
        """Calculate payoff for Asian Put option."""
        S_avg = self.S[:, 1:].mean(axis=1)
        return np.maximum(self.K - S_avg, 0)
    
    def lookback_call_payoff(self) -> np.ndarray:
        """Calculate payoff for Lookback Call option."""
        S_max = self.S.max(axis=1)
        return np.maximum(S_max - self.K, 0)
    
    def lookback_put_payoff(self) -> np.ndarray:
        """Calculate payoff for Lookback Put option."""
        S_min = self.S.min(axis=1)
        return np.maximum(self.K - S_min, 0)
    
    def price_option(self, payoffs: np.ndarray) -> float:
        """Price the option by discounting the average payoff."""
        return self.discount_factor * np.mean(payoffs)
    
    def price_european_call_mc(self) -> float:
        """Price European Call option using Monte Carlo."""
        payoffs = self.european_call_payoff()
        price = self.price_option(payoffs)
        return price
    
    def price_asian_call_mc(self) -> float:
        """Price Asian Call option using Monte Carlo."""
        payoffs = self.asian_call_payoff()
        price = self.price_option(payoffs)
        return price
    
    def price_asian_put_mc(self) -> float:
        """Price Asian Put option using Monte Carlo."""
        payoffs = self.asian_put_payoff()
        price = self.price_option(payoffs)
        return price
    
    def price_lookback_call_mc(self) -> float:
        """Price Lookback Call option using Monte Carlo."""
        payoffs = self.lookback_call_payoff()
        price = self.price_option(payoffs)
        return price
    
    def price_lookback_put_mc(self) -> float:
        """Price Lookback Put option using Monte Carlo."""
        payoffs = self.lookback_put_payoff()
        price = self.price_option(payoffs)
        return price
    
    def price_option_control_variate(self, target_payoff: np.ndarray, control_payoff: np.ndarray) -> float:
        """
        Price the option using Control Variates variance reduction technique.
        
        Parameters:
        target_payoff : ndarray
            Payoffs of the target option.
        control_payoff : ndarray
            Payoffs of the control variate (e.g., European Call).
        
        Returns:
        price : float
            Adjusted option price.
        """
        # Analytical expectation of control variate
        expected_Y = self.european_call_price_bs()
        
        # Compute means
        mean_X = np.mean(target_payoff)
        mean_Y = np.mean(control_payoff)
        
        # Compute covariance and variance
        cov_XY = np.cov(target_payoff, control_payoff)[0, 1]
        var_Y = np.var(control_payoff)
        
        # Optimal beta
        beta = cov_XY / var_Y
        
        # Adjusted payoffs
        adjusted_payoff = target_payoff + beta * (expected_Y - control_payoff)
        
        # Discounted price
        price = self.discount_factor * np.mean(adjusted_payoff)
        return price
    
    def price_asian_call_control_variate_mc(self) -> float:
        """Price Asian Call option using Control Variates."""
        european_call = self.european_call_payoff()
        asian_call = self.asian_call_payoff()
        price = self.price_option_control_variate(asian_call, european_call)
        return price
    
    def sensitivity_analysis_volatility(self, sigma_values: List[float], method: str = 'Euler-Maruyama') -> pd.DataFrame:
        """
        Perform sensitivity analysis of option prices with respect to volatility.
        
        Parameters:
        sigma_values : list or ndarray
            Array of volatility values to test.
        method : str
            Simulation method to use.
        
        Returns:
        DataFrame containing option prices for each volatility.
        """
        prices_asian_call = []
        prices_asian_put = []
        prices_lookback_call = []
        prices_lookback_put = []
        
        for sigma in tqdm(sigma_values, desc="Volatility Sensitivity"):
            # Update sigma
            original_sigma = self.sigma
            self.sigma = sigma
            # Simulate paths
            self.simulate_paths(method=method)
            # Price options
            asian_call = self.price_asian_call_mc()
            asian_put = self.price_asian_put_mc()
            lookback_call = self.price_lookback_call_mc()
            lookback_put = self.price_lookback_put_mc()
            
            # Append prices
            prices_asian_call.append(asian_call)
            prices_asian_put.append(asian_put)
            prices_lookback_call.append(lookback_call)
            prices_lookback_put.append(lookback_put)
            
            # Restore original sigma
            self.sigma = original_sigma
        
        df = pd.DataFrame({
            'Volatility': sigma_values,
            'Asian Call': prices_asian_call,
            'Asian Put': prices_asian_put,
            'Lookback Call': lookback_call,
            'Lookback Put': lookback_put
        })
        return df
    
    def plot_paths(self, num_paths: int = 10, method: str = 'Euler-Maruyama'):
        """
        Plot a specified number of simulated stock price paths.
        
        Parameters:
        num_paths : int
            Number of paths to plot.
        method : str
            Simulation method to use before plotting.
        """
        # Simulate paths
        self.simulate_paths(method=method)
        
        # Ensure num_paths does not exceed the number of simulations
        num_paths = min(num_paths, self.N)
        
        # Select the first 'num_paths' paths for plotting
        plt.figure(figsize=(12, 6))
        for i in range(num_paths):
            plt.plot(self.S[i], lw=1, alpha=0.6)
        plt.title(f'Simulated Stock Price Paths using {method}')
        plt.xlabel('Time Steps')
        plt.ylabel('Stock Price')
        plt.grid(True)
        plt.show()

# Initialize the Monte Carlo Option Pricing Engine
mc_engine = MonteCarloOptionPricing(
    S0=100,      # Initial stock price
    K=100,       # Strike price
    r=0.05,      # Risk-free interest rate (5%)
    sigma=0.2,   # Volatility (20%)
    T=1,         # Time to maturity in years
    N=1000000,     # Number of simulation paths
    ts=252,      # Number of timesteps (daily)
    seed=2024    # Random seed for reproducibility
)
# Plot the first 10 simulated stock price paths
mc_engine.plot_paths(num_paths=10000, method='Euler-Maruyama')

# Plot the first 10 simulated stock price paths
mc_engine.plot_paths(num_paths=100000, method='Euler-Maruyama')

# Initialize the Monte Carlo Option Pricing Engine
mc_engine = MonteCarloOptionPricing(
    S0=100,      # Initial stock price
    K=100,       # Strike price
    r=0.05,      # Risk-free interest rate (5%)
    sigma=0.2,   # Volatility (20%)
    T=1,         # Time to maturity in years
    N=10000,     # Number of simulation paths
    ts=252,      # Number of timesteps (daily)
    seed=2024    # Random seed for reproducibility
)

# Plot the first 10 simulated stock price paths
mc_engine.plot_paths(num_paths=100, method='Euler-Maruyama')

# Plot stock price paths using Antithetic Variates
mc_engine.plot_paths(num_paths=100, method='Antithetic Variates')


# Plot stock price paths using Quasi-Monte Carlo (Sobol)
mc_engine.plot_paths(num_paths=100, method='Quasi-Monte Carlo')

# Plot stock price paths using Milstein’s Method
mc_engine.plot_paths(num_paths=100, method='Milstein')

# Initialize the Monte Carlo Option Pricing Engine
mc_engine = MonteCarloOptionPricing(
    S0=100,      # Initial stock price
    K=100,       # Strike price
    r=0.05,      # Risk-free interest rate (5%)
    sigma=0.2,   # Volatility (20%)
    T=1,         # Time to maturity in years
    N=10000,     # Number of simulation paths
    ts=252,      # Number of timesteps (daily)
    seed=2024    # Random seed for reproducibility
)
# Plot the first 10 simulated stock price paths
mc_engine.plot_paths(num_paths=10000, method='Euler-Maruyama')

# Display Base Case Results
print("=== Base Case: Euler-Maruyama ===")
print(f"Asian Call Option Price: {price_asian_call:.2f}")
print(f"Asian Put Option Price: {price_asian_put:.2f}")
print(f"Lookback Call Option Price: {price_lookback_call:.2f}")
print(f"Lookback Put Option Price: {price_lookback_put:.2f}")

=== Base Case: Euler-Maruyama ===
Asian Call Option Price: 5.74
Asian Put Option Price: 3.33
Lookback Call Option Price: 18.31
Lookback Put Option Price: 11.75

# Control Variates for Asian Call
price_asian_call_cv = mc_engine.price_asian_call_control_variate_mc()
print("\n=== Control Variates ===")
print(f"Asian Call Option Price with Control Variates: {price_asian_call_cv:.2f}")

# Simulate using Antithetic Variates
mc_engine.simulate_paths(method='Antithetic Variates')

# Price Options with Antithetic Variates
price_european_call_antithetic = mc_engine.price_european_call_mc()
price_asian_call_antithetic = mc_engine.price_asian_call_mc()
price_asian_put_antithetic = mc_engine.price_asian_put_mc()
price_lookback_call_antithetic = mc_engine.price_lookback_call_mc()
price_lookback_put_antithetic = mc_engine.price_lookback_put_mc()

print("\n=== Antithetic Variates ===")
print(f"European Call Option Price (MC): {price_european_call_antithetic:.2f}")
print(f"Asian Call Option Price (MC): {price_asian_call_antithetic:.2f}")
print(f"Asian Put Option Price (MC): {price_asian_put_antithetic:.2f}")
print(f"Lookback Call Option Price (MC): {price_lookback_call_antithetic:.2f}")
print(f"Lookback Put Option Price (MC): {price_lookback_put_antithetic:.2f}")

# Simulate using Quasi-Monte Carlo
mc_engine.simulate_paths(method='Quasi-Monte Carlo')

# Price Options with Quasi-Monte Carlo
price_european_call_qmc = mc_engine.price_european_call_mc()
price_asian_call_qmc = mc_engine.price_asian_call_mc()
price_asian_put_qmc = mc_engine.price_asian_put_mc()
price_lookback_call_qmc = mc_engine.price_lookback_call_mc()
price_lookback_put_qmc = mc_engine.price_lookback_put_mc()

print("\n=== Quasi-Monte Carlo (Sobol) ===")
print(f"European Call Option Price (MC): {price_european_call_qmc:.2f}")
print(f"Asian Call Option Price (MC): {price_asian_call_qmc:.2f}")
print(f"Asian Put Option Price (MC): {price_asian_put_qmc:.2f}")
print(f"Lookback Call Option Price (MC): {price_lookback_call_qmc:.2f}")
print(f"Lookback Put Option Price (MC): {price_lookback_put_qmc:.2f}")

# Simulate using Milstein's Method
mc_engine.simulate_paths(method='Milstein')

# Price Options with Milstein's Method
price_european_call_milstein = mc_engine.price_european_call_mc()
price_asian_call_milstein = mc_engine.price_asian_call_mc()
price_asian_put_milstein = mc_engine.price_asian_put_mc()
price_lookback_call_milstein = mc_engine.price_lookback_call_mc()
price_lookback_put_milstein = mc_engine.price_lookback_put_mc()

print("\n=== Milstein’s Method ===")
print(f"European Call Option Price (MC): {price_european_call_milstein:.2f}")
print(f"Asian Call Option Price (MC): {price_asian_call_milstein:.2f}")
print(f"Asian Put Option Price (MC): {price_asian_put_milstein:.2f}")
print(f"Lookback Call Option Price (MC): {price_lookback_call_milstein:.2f}")
print(f"Lookback Put Option Price (MC): {price_lookback_put_milstein:.2f}")

=== Control Variates ===
Asian Call Option Price with Control Variates: 5.56

=== Antithetic Variates ===
European Call Option Price (MC): 10.47
Asian Call Option Price (MC): 5.79
Asian Put Option Price (MC): 3.35
Lookback Call Option Price (MC): 18.43
Lookback Put Option Price (MC): 11.81

=== Quasi-Monte Carlo (Sobol) ===
European Call Option Price (MC): 10.47
Asian Call Option Price (MC): 5.78
Asian Put Option Price (MC): 3.35
Lookback Call Option Price (MC): 18.34
Lookback Put Option Price (MC): 11.76

=== Milstein’s Method ===
European Call Option Price (MC): 10.60
Asian Call Option Price (MC): 5.80
Asian Put Option Price (MC): 3.42
Lookback Call Option Price (MC): 18.37
Lookback Put Option Price (MC): 11.80

# Advanced Methods Comparison
# Create a DataFrame to compare the methods
methods = ['Euler-Maruyama', 'Antithetic Variates', 'Control Variates', 'Quasi-Monte Carlo', 'Milstein']
asian_call_prices = [
    price_asian_call,
    price_asian_call_antithetic,
    price_asian_call_cv,
    price_asian_call_qmc,
    price_asian_call_milstein
]
asian_put_prices = [
    price_asian_put,
    price_asian_put_antithetic,
    np.nan,  # Control Variates implemented only for Asian Call
    price_asian_put_qmc,
    price_asian_put_milstein
]
lookback_call_prices = [
    price_lookback_call,
    price_lookback_call_antithetic,
    np.nan,
    price_lookback_call_qmc,
    price_lookback_call_milstein
]
lookback_put_prices = [
    price_lookback_put,
    price_lookback_put_antithetic,
    np.nan,
    price_lookback_put_qmc,
    price_lookback_put_milstein
]

comparison_df = pd.DataFrame({
    'Method': methods,
    'Asian Call': asian_call_prices,
    'Asian Put': asian_put_prices,
    'Lookback Call': lookback_call_prices,
    'Lookback Put': lookback_put_prices
})

comparison_df

# Sensitivity Analysis on Volatility
sigma_values = np.linspace(0.1, 0.5, 10)  # From 10% to 50%
df_sensitivity_euler = mc_engine.sensitivity_analysis_volatility(sigma_values, method='Euler-Maruyama')
df_sensitivity_antithetic = mc_engine.sensitivity_analysis_volatility(sigma_values, method='Antithetic Variates')
df_sensitivity_Quasi =  mc_engine.sensitivity_analysis_volatility(sigma_values, method='Quasi-Monte Carlo')
df_sensitivity_Milstein =  mc_engine.sensitivity_analysis_volatility(sigma_values, method='Milstein')

# Plot Sensitivity Analysis for Euler-Maruyama and Antithetic Variates
plt.figure(figsize=(12,6))
plt.plot(df_sensitivity_euler['Volatility'], df_sensitivity_euler['Asian Call'], marker='o', label='Euler-Maruyama')
plt.plot(df_sensitivity_antithetic['Volatility'], df_sensitivity_antithetic['Asian Call'], marker='s', label='Antithetic Variates')
plt.plot(df_sensitivity_euler['Volatility'], df_sensitivity_Quasi['Asian Call'], marker='o', label='Quasi')
plt.plot(df_sensitivity_antithetic['Volatility'], df_sensitivity_Milstein ['Asian Call'], marker='s', label='Milstein')
plt.xlabel('Volatility (σ)')
plt.ylabel('Asian Call Option Price')
plt.title('Sensitivity of Asian Call Option Price to Volatility')
plt.legend()
plt.grid(True)
plt.show()

Volatility Sensitivity: 100%|██████████| 10/10 [00:00<00:00, 15.12it/s]
Volatility Sensitivity: 100%|██████████| 10/10 [00:00<00:00, 19.89it/s]
Volatility Sensitivity: 100%|██████████| 10/10 [00:01<00:00,  7.51it/s]
Volatility Sensitivity: 100%|██████████| 10/10 [00:00<00:00, 12.87it/s]

# Validation: Compare MC European Call Price with Black-Scholes
# Reset simulation to Euler-Maruyama
mc_engine.simulate_paths(method='Euler-Maruyama')
price_european_call_mc = mc_engine.price_european_call_mc()
price_european_call_bs = mc_engine.european_call_price_bs()

print("\n=== European Call Option Validation ===")
print(f"Monte Carlo Price (European Call): {price_european_call_mc:.2f}")
print(f"Analytical Price (Black-Scholes): {price_european_call_bs:.2f}")
print(f"Difference: {abs(price_european_call_mc - price_european_call_bs):.4f}")

print()
# Price Options with Quasi-Monte Carlo
print(f"Monte Carlo with Quasi-Monte Carlo  Price (European Call): {price_european_call_qmc:.2f}")
# Price Options with Milstein's Method
print(f"Monte Carlowith Milstein's Method Price (European Call): {price_european_call_milstein:.2f}")

=== European Call Option Validation ===
Monte Carlo Price (European Call): 10.58
Analytical Price (Black-Scholes): 10.45
Difference: 0.1294

Monte Carlo with Quasi-Monte Carlo  Price (European Call): 10.47
Monte Carlowith Milstein's Method Price (European Call): 10.60

	Method	Asian Call	Asian Put	Lookback Call	Lookback Put
0	Euler-Maruyama	5.742512	3.331446	18.306427	11.752423
1	Antithetic Variates	5.735762	3.331488	18.251201	11.710799
2	Control Variates	5.505466	NaN	NaN	NaN
3	Quasi-Monte Carlo	5.782525	3.357352	18.318539	11.740299
4	Milstein	5.925690	3.228480	18.575486	11.521790

Monte Carlo Simulation for Pricing Asian and Lookback Options with Advanced Numerical Methods¶

Table of Contents¶

1Introduction¶

1.1. Overview¶

1.2. Objective¶

2Theoretical-Background¶

2.1 Risk-Neutral Valuation¶

2.2 Stock Price Dynamics¶

2.3 Euler-Maruyama Scheme¶

2.4 Advanced Mathematical Methods¶

3_Numerical Methods¶

Monte Carlo Simulation Steps¶

Incorporating Advanced Methods¶

Results¶

Discussion¶

Effectiveness of Advanced Methods¶

Computational Efficiency¶

Challenges Encountered¶

Insights Gained¶

Conclusion¶

Key Takeaways:¶

Future Work:¶

References¶